UNCLASSIFIED 


_ AD  NUMBER _ 

AD823871 

LIMITATION  CHANGES 
TO: 

Approved  for  public  release;  distribution  is 
unlimited.  Document  partially  illegible. 


FROM: 

Distribution  authorized  to  U.S.  Gov't,  agencies 
and  their  contractors;  Critical  Technology;  OCT 
1967.  Other  requests  shall  be  referred  to  Air 
Force  Technical  Application  Center ,  VELA 
Seismological  Center ,  Washington ,  DC.  Document 
partially  illegible .  This  document  contains 
export-controlled  technical  data. 


_ AUTHORITY 

aftac  ltr,  11  feb  1960 


THIS  PAGE  IS  UNCLASSIFIED 


168 

Revised 


FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


17  October  1967 


Prepared  For 


AIR  FORCE  TECHNICAL  APPLICATIONS  CENTER 
Washington,  0.  C. 


By 

Douglas  W.  McCowan 
TELEDYNE,  INC. 


Under 


Project  VELA  UNIFORM 


Sponsored  By 


ADVANCED  RESEARCH  PROJECTS  AGENCY 
Nuclear  Test  Detection  Office 
ARPA  Order  No.  624 


FINITE  FOURIER  TRANSFOPM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


SEISMIC  DATA  LABORATORY  REPORT  NO.  168  (REVISED) 

VELA  T/6702 

Seismic  Data  Laboratory 


AFTAC  Project  No.: 
Project  Title: 

ARPA  Order  No . : 

ARPA  Program  Code  No . : 

Name  of  Contractor: 

Contract  No.: 

Date  of  Contract: 

Amount  of  Contract: 
Contract  Expiration  Date: 
Project  Manager: 

P.  0.  Box  334, 


624 

5810 

TELEDYNE ,  INC. 

F  33657-67- C- 1313 

2  March  1967 

$  1,736,617 

1  March  1968 

William  C.  Dean 
(703)  836-7644 

Alexandria , Virginia 


AVAILABILITY 

This  document  is  subject  to  special  export  controls  and  each 
transmittal  to  foreign  governments  or  foreign  national  may  be 
made  only  with  prior  approval  of  Chief,  AFTAC. 


V 


This  research  was  supported  by  the  Advanced  Research 
Projects  Agency,  Nuclear  Test  Detection  Office,  under  Project 
VELA-UNIFORM  and  accomplished  under  the  technical  direction 
of  the  Air  Force  Technical  Applications  Center  under  Contract 
F  33657-67-C-1313. 

Neither  the  Advanced  Research  Projects  Agency  nor  the 
Air  Force  Technical  Applications  Center  will  be  responsible 
for  information  contained  herein  which  may  have  been  supplied 
by  other  organizations  or  contractors,  and  this  document  is 
subject  to  later  revision  as  may  be  necessary. 


TABLE  OF  CONTENTS 


Page  No. 


ABSTRACT 

1.  INTRODUCTION  1 

2.  THE  FINITE  AND  DISCRETE  FOURIER  TRANSFORMS  1 

3.  TWO-AND  THREE-DIMENSIONAL  FOURIER  TRANSFORMS  5 

4.  ALGEBRAIC  DISCUSSION  7 

5.  HIGH-SPEED  CORRELATIONS  AND  CONVOLUTIONS  12 

REFERENCES 

APPENDIX  A  -  PROCEDURES 
APPENDIX  B  -  PROGRAM  LISTINGS 


APPENDIX  C  -  PROGRAM  WRITE-UPS 


ABSTRACT 


The  theory  of  finite  Fourier  transforms  is  developed 
from  the  definitions  of  infinite  transforms  and  applied  to 
the  computation  of  convolutions,  correlations,  and  power 
spectra.  Detailed  procedures  for  these  computations  are 
given,  including  listings  and  writeups  of  FORTRAN  subroutines. 
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INTRODUCTION 


I 


A 


f 


For  the  past  several  months,  E.  A.  Flinn,  j.  F.  Claerbout, 
and  I  have  been  examining  some  practical  and  computational 
aspects  of  the  theory  of  Fourier  transforms.  These  efforts  have 
resulted  in  a  set  of  programs  for  performing  operations  on  time 
series  based  on  the  Cooley-Tukey  (References  1,2)  hyper-rapid 
Fourier  transform  method.  Using  this  method,  computations  on 
seismic  array  data  such  as  the  calculation  of  convolutions, 
correlations,  spectra,  and  digital  filters  have  been  speeded  up 
by  factors  of  three  or  four  and  sometimes  even  ten.  The  purpose 
of  this  report  is  to  communicate  these  results  i>  a  straight¬ 
forward  manner  and  to  offer  some  motivation  for  their  derivation 
as  well  as  for  future  efforts  in  this  rrea.  Writeups  and 

listings  of  the  programs  discussed  here  are  included  as  appendi¬ 
ces  to  this  report. 


2*  THE  FINITE  AND  DISCRETE  FOURIER  TRANSFORMS 

In  the  case  of  continuous  data  of  infinite  length,  the 
Fourier  transform  pair  is  usually  written  as: 


A(u>)  =  (2tt)  2  J  f(t)e'iu,tdt 


f(t)  =  (2tt)  2  J  A(ui)e3 


The  first  of  these,  going  from  time  to  frequency,  is  referred 

I 
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to  as  the  direct  transform  and  the  other  as  the  inverse  trans¬ 
form.  Sometimes  the  direct  transform  is  written  with  a  factor 
of  1  in  front  of  the  integral  and  the  inverse  with  a  factor  of 

I/2 n  ’  These  are'  of  course,  equivalent  to  the  above  definition 
Quantities  of  interest,  such  as  spectra,  etc.,  involve  magni¬ 
tudes  or  squares  of  one  transform  and  the  factor  must  be 
inserted  or  taken  out,  depending  on  which  definition  is  used, 
to  preserve  the  proper  dimensions. 

Two  drawbacks  of  these  definitions  for  digital  compu¬ 
tations  are  apparent:  First,  the  integrals  must  be  approxi¬ 
mated  by  sums  in  the  digital  computer ,  which  implies  that  both 
transforms  involve  sampled  variables.  Second,  the  infinite 
limits  on  the  sums  are  impossible  in  practice.  Clearly,  these 
oums  must  be  truncated,  as  they  do  not  in  general  converge  over 
a  finite  interval.  As  a  result  Fourier  transforms  as  such  are 
never  really  computed  by  a  digital  computer.  Instead,  the 
complex  samples  of  a  direct  transform  are  approximated  by  the 
cosine  and  sine  coefficients  of  Fourier  series  representation 
of  the  input  data.  The  definitions  for  these  are: 

00 

if  X(t)  =  I  [an  cos  <"nt/T)  +  bn  sin  (nnt/T)]  ,  (2) 


then  aQ  =  i  J  x(t) 
o 


b  =  0 
o 


2  r 

an  “  T  J  x<t)  cos  (nnt/T)  dt 


2 


sin  (rnit/T)  dt 


bn  *  f  J  x(t> 

o 

If  N  samples  of  the  data  are  taken  at  equally  spaced 
intervals  At  =  T/N ,  the  integrals  (3)  becomes  sums  and  the 
frequency  sum  in  (2)  goes  from  DC  to  the  folding  frequency, 
i.e.,  k  =  0  to  N/2T  .  The  equations  are  then  written  as; 


x(j)  -  £  [ak  cos  (  2tt  jk/N )  +  bk  sin  (2njk/N)J 


3o  "  N  ZX(3) 


b  =  0 
o 


3k  "  N  I  X(j)  cos  ( 2n jk/N ) 
j“0 


bk  ~  n  I  x(j)  sin  ( 2n jk/N ) ,  (5, 


where  t  has  been  replaced  by  jAt  .  By  now  def 


ining: 


A(K)  =  2  (ak  '  1  V  '  A(0>  =  a0  (6) 

and  realizing  th  it  a  real  time  series  contains  only  real  points,  we 
can  write  (4)  ass 


x(j)  s  Yj  exP  (2rri  jk/N) 


A  groat  deal  of  symmetry  between  the  two  transforms  can  be 
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preserved  if  the  sum  in  (7)  is  summed  up  to  N-l  .  Redundant 
points  in  the  spectrum  are  included  (since  the  transforms  are 
periodic)  but  the  computational  procedures  are  simplified. 

It  is  also  convenient  to  split  the  factor  of  1/N  appearing  in 
(5)  into  two  factors  of  1//N  ,  one  in  front  of  each  transform. 
By  defining  a  complex  number? 

w  =  exp  (2rri/N)  ,  (8) 

the  two  transforms  can  now  be  written  as 

N-l 

MM  *  7f  I  £<j)  w'ik  <«> 

j-0 

N-l 

f(j)  =  7f  ^  a(m  wjk  uo) 

k=0 

It  can  be  shown  that  the  set  of  direct  Fourier  trans¬ 
form  points,  between  DC  and  the  folding  frequency,  contains 
the  same  amount  of  information  as  the  real  data  series:  The 
transform  includes  N/2  distinct  points,  which  with  the  DC  term 
makes  a  total  of  N/2  +  1  complex  points.  Equation  (9)  shows 
that  both  the  DC  and  the  folding  frequency  point  are  purely 
real?  thus,  the  Fourier  transform  contains  (N/2-l)*2+2*l  num¬ 
bers.  This  is  exactly  the  same  amount  of  information  con¬ 
tained  in  the  real  time  series.  It  also  suggests  that  the 
existence  of  one  transform  should  imply  the  existence  of  the 
other . 

If  there  are  N/2+1  non-redundant  points  in  the  direct 
transform,  then  the  sampling  interval  in  frequency  must  be 
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(N/2T )/ (N/2 )  -  1/T  .  Thus,  the  product  of  the  time  and  frequency 


variables  is; 


iuut  =  i(2nk"— )  (j”)  =  jk 
T  JN  N  J 


This  equation  relates  the  arguments  in  the  two  exponentials, 
one  in  the  continuous  transform  and  the  other  in  the  finite 
transform  (Equations  1,  9,  and  10). 

3.  TWO-  AND  THREE-DIMENSIONAL  FOURIER  TRANSFORMS 

Two-  and  three-dimensional  direct  Fourier  transforms 
are  seen  to  be 


V1  n2-i 


(kl'k2)  /NT  Z  Z  X^l'^2*  W1  ^  1  W2  2 


12.  „  .  „ 
3r°  V° 


N  -1  N-l  N  -1 
12  3 


A(k1#k2,k3)  = 


N  N  Z  Z  Z  X ^1 '  ^2 '  ^ 3^  W1  ^  1 


1  2  3  VO  j2=o  j3=o 


w2  j2k2  w3_:,3k3  (13) 


We  can  break  up  Equation  (12)  as  follows: 


n2-i 


A(ki,k2)-/_  £  B(kl'j2)  W2  j 2  2 


V° 


5 


This  calculation  requires  one-dimensional  transforms;  we 
have  defined 


V1 


B(k 


i'j2»  't sr  l  x<w  wi'3lkl 


1  v° 


(15) 


which  requires  N  one-dimensional  transforms.  Thus,  N  +  N 

2.  12 

one-dimensional  transforms  are  required  to  compute  the  single 

two-dimensional  transform. 

We  can  break  up  Equation  (13)  as  follows: 


«3-l 


a(k1.k2,k3)  yjj---  £  C(k1,k2,.i3)  v#3  :i3k3 


(16) 


3  v° 


which  requires  N^N.^  one-dimensional  transforms;  we  have  de¬ 
fined 


V1 


^ (^^ '^2 ' ^  3^ 


1  2 


jj*0 


V1 


/nTnT  I  I 


x(j1»j2#j3)  wx  ^ 1^1  w2  ^2k2  (11) 


“1  =■ 
J2 


which  requires  two-dimensional  transforms.  Thus,  one¬ 

dimensional  transforms  and  N,  two-dimensional  transforms  are 

j 

needed  to  compute  the  single  three-dimensional  transform. 


4 .  ALGEBRAIC  DISCUSSION 

Equations  (9)  and  (10)  suggest  a  more  elegant  and  compact 
way  to  write  the  two  transfor  s.  We  define  the  vector  A  as 

the  transform  with  elements  (A)^  =  A(k),  and  define  the  vector 

F  as  the  time  series  with  elements  (F)^*F(j)  .  The  process 

of  transforming  is  seen  to  be  equivalent  to  matrix  multipli- 

cation  by  a  matrix  W  whose  elements  are  (W)  =  wJ 

jK 

A  =  W+  f  (18) 


arid  f  «  WA  ,  (19) 

where  the  dagger  indicates  Hermitian  conjugation.  Substituting 
(19)  into  (18)  gives  the  following  important  identity: 

WW+  *  W+W  =  I  .  (20) 

This  is  the  definition  of  unitarity  for  the  transformation  W 
It  is  a  generalization  of  orthogonality  for  complex  matrices 
and  assures  Parseval ' s  theorem: 

A+A  =  f+f  .  (21) 

W  preserves  "length"  between  the  two  domains.  The  identity  is 
actually  proved  by  writing  out  the  terms  in  the  product: 


^  £  [exp(2ni/N)  J  j^exp(-2ni/N)  J 


=  6: 


7 


or 


< 


N-l 

I  £  W-<H<> 

m=0 


(22) 


This  last  important  relation  is  seen  to  be  true  by  the  use  of 
a  phase  diagram; 


The  Cooley-Tiikey  method  factors  the  W  matrix,  if  its 
order  is  a  power  of  two,  into  L  +  1  sparse  matrices,  where  L 
is  the  power  of  two: 


L-l 


S 


0 


Multiplying  L  +  1  times  by  these  sparse  matrices  can  in  some 
cases  reduce  the  computing  time  by  many  tens  of  times.  This 
factorization  is  proved  by  Good  (3)  and  organized  for  compu¬ 
tation  by  Rader  (4). 

For  the  case  N  —  8  the  W  matrix  is: 


8 


1_  1 


1 

1 

l 

1 

1 

1 

1 

1 

w 

2 

w 

3 

w 

4 

w 

5 

w 

6 

w 

7 

w 

2 

w 

4 

w 

6 

w 

8 

w 

10 

w 

12 

w 

14 

w 

3 

w 

6 

w 

9 

w 

12 

w 

15 

w 

18 

w 

21 

w 

4 

w 

8 

w 

12 

w 

16 

w 

20 

w 

24 

w 

28 

w 

5 

w 

10 

w 

15 

w 

20 

w 

25 

w 

30 

w 

35 

w 

6 

w 

12 

w 

18 

w 

24 

w 

30 

w 

36 

w 

42 

w 

7 

w 

14 

w 

21 

w 

28 

w 

35 

w 

42 

w 

49 

W  | 

(22.1) 


The  L  +  1  =  4  transformations  are  graphically  illustrated 
by  the  following  diagram  in  Rader's  notation  (Reference  4): 


INDEX  BINARY  ARRAY 


x(o)r - - ►A  ( 0) 

x  ( i  H  V  •  f® Ml) 

X(2)r^'/f®AfeT>©--VHk»A(2) 

X(3)ry'/^o)Z_^)z^e)v  /\  A( 3) 

\F.  r  l  2  —  v —  v  .  \J 


A  (4  ) 


X(4,4^77^@v“k^'^0'^  /^A(4) 

X  (  5  )4y^0  *  -  a  (  5 ) 

X(6  \(6) 

X(7)^ )4d  ►  A  ( 7 ) 


X(6> 

X(7) 


A  ( 6) 
A  ( 7 ) 


RE¬ 
VERSED  INDEX 


Ec  :h  solid  line  represents  a  multiplication  by  w  to  the 
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power  indicated  in  the  circle,  and  each  dotted  line  represents 
a  simple  addition  into  that  element  of  the  array..  No  ad¬ 
ditional  storage  is  used  by  this  process.  The  results  of  each 
transformation  are  stored  on  top  of  the  original  data,  and 
the  last  transformation,  which  is  a  simple  interchange,  gives 
the  desired  Fourier  transform.  Note  also  that  the  succession 
of  numbers  in  the  circles  is  the  bit-reversed  representation 
of  the  sequence  of  indices  in  order.  They  can  be  stored  in  a 
table  or  generated  successively  by  a  reverse-add  procedure. 

For  reasons  of  space  and  simplicity,  we  have  chosen  the  latter 
route . 

The  S  matrices  are: 


f1 

0 

0 

0 

1 

0 

0 

1. 


0  0  0 
10  0 
0  10 
0  0  1 


w  0  0  0 

0  w^  0  0 

0  0  w^  0 


0 


0  0 


0 

w 


0  0  0 
10  0 
0  10 
0  0  1 
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V 


1 

rr 
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HIGH-SPEED  CORRELATIONS  AND  CONVOLUTIONS, 


By  computing  Fourier  transforms  with  this  finite  Fourier 
series-like  method  an  important  condition  is  put  on  the  time 
series.  As  in  regular  Fourier  series  the  input  is  assumed  to 
be  periodic  with  period  T  and  the  integrals  or  sums  are  com¬ 
puted  over  a  single  period.  There  is  also  the  effect  of 
cutting  off  the  spectrum  at  the  folding  frequency.  Sines  and 
cosines  of  finite  wavelength  will  repeat  again  outside  the 
region  of  interest.  This  fact  in  itself  is  not  bothersome  but 
becomes  a  serious  complication  in  the  computation  of  convo¬ 
lutions  and  correlations.  Convolutions  and  correlations  as 
usually  computed  assume  the  time  series  to  be  zero  outside  the 
region  of  interest.  Therefore,  the  integrals  or  sums  in  com¬ 
puting  them  are  summed  out  only  over  the  non-zero  terms.  When 
multiplying  together  two  finite  Fourier  transforms  (or  the 
complex  conjugate  of  one  times  the  other)  the  periodicity  of 
the  time  series  means  that  elements  which  have  been  shifted 
past  the  end  of  a  period  reappear  at  the  beginning.  This 
process  is  called  circular  convolution  or  correlation  and  its 
effects  are  unavoidable  when  straightforwardly  computing  lagged 
products  with  finite  Fourier  transforms.  This  is  illustrated 

below; 

X1  =  (3'  °'  ~1'  2) 

X2  =  (  -2,  2,  -1,  3  ) 

Rc  =  (i,  -l,  3,  5)  for  100%  positive  lags; 

12  * 

-  (1,  5,  3,  -1)  for  100%  negative  lags. 
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Circular  convolution  is  therefore  written? 


whore  the  upper  limit  on  the  sum  simulates  the  desired  zeros 
in  the  time  series  outside  the  region  of  interest.  This  is 
illustrated  below: 

Xx  =  (3,  0,  -1,  2) 

X2  =  (-2,  2,  -1,  3) 

R12  =  3'  ~3'  f°r  100^  positive  lags; 

=  (1,  -4,  6,  -4)  for  100%  negative  lags. 

T 

The  finite  Fourier  transform  of  this  R  is  thus  not  the 
product  of  the  two  individual  transforms.  However,  by  filling 
zeros  into  the  second  half  of  each  data  series  nd  computing 
their  transforms  out  to  twice  their  actual  length,  a  good  esti 
mate  of  the  spectrum  may  be  obtained.  In  addition,  the 
negative  lags  in  the  correlation  appear,  thus  giving  a  more 
mathematically  satisfying  result.  This  is  illustrated  below: 

Xl  =  (3,  0,  -1,  2,  0,  0,  0,  0) 

X2  =  (-2,  2,  -1,  3,  0,  0,  0,  0) 

R?„  =  (1,  3,  -3,  9,  0,  -4,  6,  -4)  for  100%  positive  lags 


The  two  modified  transforms  thus  are: 

2T-1 

F. (k)  =  Y  x. (t)  w"tk  x. (t)  =  0,  T  <  t  <  2T-1 

l  o  l  l 

t=0 

2T-1  2T-1 

S..(k)  =  F, (k)*  F.(k)  =  V  x.(t)  wtK  Y  x  (t)  w  T 
lj  1  j  l  <->  j 

t=0  r=0 
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2T-1 


'ij(s)  *  I  ri(k>*  FjOO  « 


2T-1  2T-1 


2T-1 


I  I  Xi(t>  x-j(T)  £  wk(t+S"T) 


t=0  T=0 


Now  from  (22)  the  last  sum  becomes  a  Kronecker  delta  function 
and  the  other  sum  is  collapsed  to  give: 


2T-1 


"ij15’  *  I  xi(t>  *}<*+»>  -  Rj .(s) 


The  last  equality  following  from  the  original  assumption  that 
Xi*t*  °>  T  <  t  <  2T-1  .  Transient  correlations  for  100K 

lags  are  therefore  computed  by  forming  the  absolute  product  of 
two  transforms,  each  computed  out  to  twice  the  length  of  the 
original  data  series  with  zeros  filled  into  the  second  halves. 

Non-circular  or  transient  convolutions  are  computed  In 
much  the  same  way.  except  that  the  transforms  have  to  be  com¬ 
puted  out  to  a  length  equal  to  the  sum  of  the  lengths  of  the 
time  series  and  the  filter,  with  the  appropriate  number  of 

zeros  filled  into  each.  The  convolution  theorem  is  proved  in 
the  same  fashion. 


T+S-l 


l(k)  =  I  a(T)  w_Tk  a(r)  =  0,  S£t<T+S-1 


15 


x(t) 


x(t)  =  0 


T+S-l 

X-OO  -  l 

t  =  0 


T<t<T+S-l 


T+S-l 

T+S-l  T+S-l 

T+S-l 

z 

A(k)  X(k) 

ku 

w 

=  Z  Z  a,T) 

x(t)  J 

k (u-T-t) 
w 

k=0 

T  =  0  t  =0 

U 

k=0 

T+S-l 

T+S-l 

z 

A  (k )  x  (k ) 

ku 

w 

=  ^  a(t )  x (u-t )  = 

y(u) 

(27) 

k=0 

T  =  0 

Where  y(u)  is  now  the  "filtered"  output  of  the  filter  a  acting 
on  X.  Convolutions  are  therefore  computed  by  forming  the 
product  of  the  twc  transforms,  each  computed  out  to  a  length 
equal  to  their  sum  with  zeros  filled  into  the  extr.  lengths. 
Detailed  procedures  for  these  commutations  are  listed  in 
Appendix  A . 
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APPENDIX  A  -  PROCEDURES 


FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


\ 


sp 


PROCEDURE  FOR  CAI,CULAT tN(!  AN  AtrTO- SPECTRUM  AND  AN  AUTO¬ 
CORRELATION _ _ _ _ 

DIMENSION  X(2*LX+2),  CX(LX+1) 

EQUIVALENCE  (X.CX) 

TYPE  COMPLEX  CX,  CONJG 
LX  -  2**N 


l 


1)  Erase  2+LX+2  points  in  Xr  the  extra  complex  point  Is  needed 
by  COOLER  to  return  the  point  at  the  folding  frequency. 

2$  Read  the  data  channel  Into  X(l)  through  X(LX) . 

3)  CALL  COOLER ( N+ ). ,CX) .  The  Fourier  transform  of  X  and  the 
necessary  zeros  on  the  end  of  the  data  Is  now  stored  In 
CX,  LX+1  complex  points  long,  representing  frequencies 
between  DC  and  the  folding  frequency. 

4)  Go  through  the  LX+1  complex  points  In  CX,  andt 

CX(I)  -  rcONJG(CX(I))*CX(n)/LX 


that  la, 

*  Re[CX(I)  1  -  (RelCX(I)l2  +  IrnfCXd )  ]2)/LX 

ImrcX(I) ]  -  0.0 

The  auto-spectrum  Is  the  real  part  of  CX,  purely  real 
and  LX+1  points  in  length. 

5)  To  get  the  auto-correlation,  fill  In  the  other  LX-1 
complex  points  in  CX  as  required  by  COOL  for  inverse 
transforms,  and  call  COOL) 

DO  1  I  -  1 , LX-1 
1  CX(LX+ 1+1 )  -  CX (LX-I+1 ) 

CALL  C00L(N+1 ,CX, +1.0) 

The  auto-correlat Ion  is  in  the  real  part  of  CX,  purely  real 
and  2*LX  points  in  length. 


NOTE:  CX  must  be  dimensioned  2*I.X  if  the  auto-rorreldtion 
is  to  be  computed. 

-  A-l  - 


PROCEDURE  FOR  CALCULATING  A  CROSS  SPECTRUM  AND  A  CROSS -CORRELATION 


DIMENSION  X(2+LX+2 ) ,  CX(LX+1),  Y(2*LX+2),  CY(LX+1) 
EQUIVALENCE  (X,CX).  (Y,CY) 

TYPE  COMPLEX  CX.CY 
LX  -  2**N 

1)  Erase  2*LX+2  points  i.n  both  X  and  Y. 

2)  Read  channel  1  into  X  and  channel  2  into  Y. 

3)  CALL  COOLER(N+l  ,X) 

CALL  COOLER(N+l.Y) 

4)  Go  through  the  LX+1  complex  points  and  overlay  CX  (or  CY) 
with* 

CX(I)  -  [CONJG(CX(l))*CY(I)]/LX 
that  is. 

Re[CX(I) 3  -  (Re[CX( I) ]*Re (CY( I) ]+ImCcX(l) 1+Im[CY(I) ])/LX 

ImfCX(I) ]  -  (Re[CX(in*Im[CY(I)]-ImrcX(I)l*Re[CY(l)])/LX 

The  cross-spectrum  between  channel  1  and  channel  2  (Which 
is  the  complex  conjugate. of  the  cross-spectrum  between 
channel  2  and  channel  1)  is  now  in  CX,  LX+1  points  in 
length.  The  co-spectrum  is  in  the  real  part  of  CX  and  the 
quad-spectrum  is  in  the  imaginary  part  of  CX. 

5)  To  get  the  cross-correlation,  fill  in  the  other  LX-1  points 
in  CX  and  call  COOL: 

DO  1  I  “  1 , LX-1 

1  CX ( LX+ 1+ 1 )  =  C0NJG(CX(LX-I+1) ) 

CALL  C00L(N+1 ,CX ,+1.0) 

The  cross-correlation  is  in  the  real  part  of  CX,  purely 
real  and  2*LX  points  in  length. 

NOTE:  CX  must  be  dimensioned  2*LX  if  the  cross-correlation 
is  to  be  calculated. 
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PROCEDURE  FOR  CALCULATING  THE  COMVm.irrr™ 


OF  TWO  SERThS 

DIMENSION  X  ( L+  2 )  ,  CX(»jL+l),  F(L+2>,  CP(«jL+l) 

EQUIVALENCE  (x,CX).  <p,CF* 

TYPE  COMPLEX  CX.CF.CONJG 
L  “  2**N 

L  here  i.  the  next  power  of  2  larger  than  LX+LF,.the  combined 
length  of  the  data  and  the  filter. 

1)  Erase  L+2  points  in  X  and  F. 

2)  Read  the  data  into  X(l)  through  X(LX)  and  the  filter  impulse 
response  into  F ( 1 )  through  F(LF). 

3)  CALL  COOLER (N,CX) 

CALL  COOLER (N,CF) 

4)  Go  through  the  *iL+l  complex  point*  in  CX,  andt 

CX(I)  -  rcX(I)*CF(I)  (/LX 
that  is, 

Re|CX(I)1  -  (RefCX(I)  l*Re(CF(I)  l-Im|CX(I) 1*Im|CF(I) ])/LX 

ImfCX (1)1  -  (Rercx(I)  |*ImfCF(I)l+RercF(l)1*lmrcx(l)  |)/LX 

The  Fourier  transform  of  X  convolved  with  r  is  now  in  CX. 

5)  Fill  in  the  re ’t  of  the  points  in  CX  as  needed  by  COOL, 
and  transform  back.  Note  again  that  if  the  actual  convo¬ 
lution  is  desired  instead  of  the  Fourier  transform,  CX 
must  be  dimensioned  L. 

DO  )  t  *  1,  'jL-l 

1  CXC/L+I+I)  «  CONJGTcX (!jl,  -I4J)1 
CALL  COOLiN ,CX ,-1.0) 

The  convolution  of  X  with  F  i s  now  in  the  real  part  of  CX, 
purely  real,  and  !.X*LF-1  points  in  length. 


APPENDIX  B  -  PROGRAM  LISTINGS 

FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


aaoaanaaanaaaaananaaaannnnonanaannnnaannannnnnnnannanaa 


SUBROUTINE  CUOUN,X#SlGNl> 


HYPEH-RAPID  FUURIPR  TRANSFORM  USING  COOLEY-TUKEY  ALGORITHM 


SEISMIC  DATA  LABORATORY*  ALEXANDRIA*  VA,  PROGRAMMED 
26  FEBRUARY  1V66  RY  J«  F «  CLAfcRBOUT  (HIT)*  0,  W,  MCCOWAN, 

E «  A «  FLINN#  AND  JI  GIBSON  I  TELEDYNE ) 

X  IS  A  COMPLEX  ARRAY  USED  FOR  THE  DATA  SERIES  AND  THE 

N 

transform  -  the  number  of  elements  of  x  is  l  >  2  • 

sign  •  -i.o  fur  direct  iourier  transform  and  ♦1.0  for  inverse 

FOURIER  TRANSFORM  (BUT  SEE  BELOM  FOR  ARRANGEMENT  OF  DATA  FOR 
INVERSE  TRANSFORM). 

FOR  DIRECT  TRANSFORM*  ON  INPUT  THE  REAL  PART  OF  X  CONTAINS  THE 
DATA  SERIES  AND  ThE  IMAGINARY  PART  OF  X  IS  ZERO.  BN  RETURN* 
THE  FOURIER  COSINE  SERIES  EXPANSION  OF  THE  DATA  IS  IN  THE  REAL 
PART  OF  X*  AND  THE  FOURIER  SINE  SERIES  EXPANSION  IS  IN  THE 

N-l 

IMAGINARY  PAHT  OF  X.  EACH  CONTAINS  ONLY  l  *  1  NONREDUNDANT 
POINTS.  THE  COSINE  EXPANSION  IS  SYMMETRIC  ABOUT  POINT  NUMBER 

N*1 

2  *  1  AND  The  SINE  TRANSFORM  is  ANTISYMMETRIC  ABOUT 

THIS  POINT, 

FOR  EXAMPLE  "  N  «  3  AND  DATA  •  <0.*l.*0}*0.*0.*0.*|.»0v>, 

THEN  REAL  PART  OF  X  «  U,*1.,0,,0.»Q.*0; ,0.#0. >  AND  IMAGINARY 
PART  OF  X  •  (Q,*0.«0t»0«*0**0,,0,*0*)  ON  INPUT* 

ON  RETURN*  REAL  PART  OF  X  ■  U , 000* *7071, 0* *•* 7071*-1 , 000* 
•*70/1*0, .♦./Oil)  AND  IMAGINARY  part  OF  X  •  (0««*,?071, 

•1. 000*-, 7071*0.#. 7071*1. 000,,7071>.  POINT  NUMBER  1 
CORRESPONDS  TU  ZERO  FREQUENCY*  POINT  NUMBER  9  CORRESPONDS 
TO  PI*  THE  FOLDING  FREQUENCY, 

TO  DO  AN  INVERSE  TRANSFORM*  THE  COSINE  AND  SINE  SERIES  MUST  BE 

N»1 

FOLDED  OVER  ABOUT  POINT  NUMBER  2  *  1  BEFORE  CALLING 

COOL  WITH  SIGN  ■  ♦l.O*  SUBROUTINE  FTPACK  CAN  BE  USED  TO  DO 
THIS  FOR  YOU*  CONVERTING  AMPLITUDE  AND  PHASE  BACK  TO 
SINE  AND  COSINE  IF  NEED  BE* 

-N 

THERE  IS  A  SCALE  FACTOR  OF  2  WHICH  COOL  DOES  NOT  APPLY, 


THE  USER  CAN  APPLY  THE  SCALE  FACTOR  EITHER  TO  THE  DIRECT  OR  TO 

-N/2 

THE  INVERSE  TRANSFORM*  OR  APPLY  A  SCALE  FACTOR  OF  2  TO 


C  FOR  EXAMPLE#  GIVEN  the  input  uata  as  above#  the  two  STATEMENTS 

C  CALL  COOL ( 3#  X#  *1 , o  ) 

C  CALL  COOL(3»X#*1(0) 

C  WOULO  CHANGE  HfiAb  PART  OF  X  TO  < 0 # #  8 , , 0*  ,  0  . ,  0  .  # 0  .  *9 . # 0 « »  AND 

C  IMAGINARY  pant  OF  X  TO  <0*»0'#0»#0t,0,,0((0.,0*>. 

c 

c 

c 

DIMENSION  XU>#  INTU6)*fl<2> 

TYPE  COMPLEX  X#U#W,HOLD 
equivalence  <q#w> 
c 

c  INITIALIZE 

C 

LX  •  2**N 
PI2-0. 203189306 
FLX  •  LX 

FlXPI2»S1GNI*RI2/FLX 
DO  10  1 Rl#N 
10  INTI  I  )  ■  2**<N»1> 

C 

C  LOOP  OVER  N  LAYERS 

C 

DO  40  LAYER  •  1#  N 
NBLOCK  •  2**(LATfaR«l| 

lblock»lx/nblqck 

LBHALF  •  LBLOCK/2 
C 

C  START  SERIES  AND  LOOP  OVER  BLOCKS  IN  EACH  LAYER 

C 

NW  »  0 

DO  40  IPL0CK"1«NBL0CK 
LSTAKT  ■  LBLOGKM  JBLOCK’l) 

C 

C  COMPUTE  W  ■  CbXPI2»*PI*NW*SlQNI/LX) 

C 

ARG«FL0ATF«NW)*FLXPI2 
GUT  ■  COSFTARQ} 

0(21  •  SlNFiARGI 

C 

C  THIS  CAN  BE  SPEEDED  UP  BY  USING  A  TABLE  OF  COSINES 

C 

C 

C  CQMPUIE  ELEMENTS  FOR  BOTH  HALFS  OF  EACH  BLOCK 

C 

DO  20  1«1#LBHALF 

J  ■  ulstart 

K  •  U4LBHALF 
Q  ■  X(KI*W 
X  (K  )  ■  X ( J ) »Q 
X(J)  ■  X(J)*Q 
20  CONTINUE 
C 

C  BUMP  UP  SERIES  BY  TWO  (NOT  ONi) 

C 

DO  32  I ■2#N 


1 1  »  1 

LL«1NT( 1 » .ANU.NW 


C 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


THIS  LOGICAL  OPERATION  IS  A  MASK  TO  DETECT  A  ONE  If* 

THE  AHPROPRIAlfi  BIT  P09ITI0N  OF  NW,  THIS  STATEMENT  WILL  NOT 
WORK  UN  IBM  FORTRAN  STS  I fcMS  # 

IF ( LL )31«  3l»  30 

30  CONTINUE 

NW  a  NW«INT(  I  ) 

32  CONTINUE 

31  CONTINUE 

NW  »  NWMNTCIP 
40  CONTINUE 

START  SERIES  10  BEGIN  FINAL  REPLACEMENT 
NW  ■  0 

00  50  K»l#LX 

CHOOSE  CORRECT  INDEX  AND  SWITCH  ELEMENTS  IF  NOT  ALREADY 

sw itched 

NWlaNW«l 

IF (NW1«K)55»99»6U 
60  HOLDaXINWll 
X(NWIMXIK) 

X  c  K  >  ■  HOLD 
55  CONTINUE 


BUMP  Up  SERIES  BY  ONE 


00  70  I ■!# N 
II  ■  I 

LL-1NTU  )  .AND.NW 
IF (LL  ABO*  80#  70 
70  NW  ■  NWINT(I) 

80  NW  •  NW+ 1  NT (  1 1 ) 
50  CONTINUE 
RETURN 
END 


\ 

\ 
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9UHR0UTINE  COQLCU|im|NT,IOT.L»F.X> 

0|M£NSION  F(l)»X«2,l),LAU(12/> 

DIMENSION  F(N,L>*X(2*I?E$T>*LAB(127) 

MULTICHANNEL  CONVOLUTION  ROUTINE  F OH  .UPEB  DATA 

INf  |S  THE  INPUT  SUBSET  TAPE  6f  DATA  CHANNELS 
IOT  IS  TNfc  OUlRUT  SUBSET  UPE  OF  UATA  CHANNELS 
L  IS  THE  NUMBER  OF  FILTER  POINTS  F OH  EACH  CHANNEL 
F  IS  I  HE  f JLIER  MATRIX 

X  IS  A  WORKING  ARRAY  CONTAINING  AT  LEAST  2*ITEST  POINTS 
ITEST  IS  THE  NEXT  POWER  OF  I  HO  LARGER  THAN  LX*L 

C  D.M.MOCOWAN  JULY  1966 

C 

REWIND  I  NT 
REWIND  IOT 
REAUdNHLAB 
N*LAB(2) 

LX«L*B< J> 

ISUM*LX*L 
LAB(  «J)«LX-CL-1> 

writeciohlah 

DO  1  lNO"l,l«J 
I TES1 *2** IND 
II'(|S0M*1TEST)2*2j1 

2  NCOOL-INU 
00  TO  3 

1  CONTINUE 

PRINT  lOUO.LX.L 

1000  FQRMATC59H1BAD  NEWS*  ERROR  IN  COOLCON,  DATA  PLUS  FILTER  TOO  LONG  L 
IX*  «|6*5H*  l*  *161 
STOP 

3  CONTINUE 
IT02-ITEST/2 

I T02H2* I T  02*2 

DO  10  IN>1,N 

CALI  ERASE«2*ITEST#X) 

READ!  INT > (X(1,M>»m*1#LX) 

00  11  ILB1*L 

U  X(2# 1L>»F ( |N*|IL-1)*N) 

CALL  COUL(NCQOL»X*-1,OJ 
X (1«1)*X (1*1>«X(2*1)/|TEST 
X(2«l)*0«0 
DO  20  IL*2« I T02 

9AVE*(XI1* ITEST«ILf2)*X(2« I  TEST- 1 L*2 > *X ( 1# 1 L > *X ( 2* I L > >/ (2*1  TEST ) 

lJ<2,IU*(X(l»,TEST.,L*2»*#2.X(2,|TEST-a»2l**2»X(l»IL>*A2*X(2,,Ll* 
20  X(1«ILI*SAVE 

X(l« I I02*1)*X(1«1To2+1)*X(2-1T02*1)/1TEST 
X(2. 1T02*1)*0,0 
00  30  IL*IT02P2* (TEST 
X(1,1L)«X<1, ITESI.IL+2) 

30  Xt2«IL>B‘*X(2*IlEST-lL+2> 

CALL  COOL(NCU0L*x«*1,0) 


10  WRITEdOl  )(X(1»MI#m*l,lX) 
end  f  ile  i o t 

HEWIND  IUT 

rewind  INT 

HE  I  URN 
END 


oonono 


SUBROUTINE  CQQLfcH ( N»  X ) 

DIMENSION  XA2*l*H 

FOUR  I bR  TRANSFORM  OF  A  HgAL  DATA  SERIES 
N  MUSI  BE  LESS  THAN  OR  EQUAL  TO  14 

DIMENSION  y  (2,1) 

M  ■  N-l 
L  ■  2**M 

Call  coukm,x,-i.o> 
f  ■  3,14159265J/FL0ATF(U 
SAVE-XU) 
mixxaHxu) 

X(l,Ln>*SAVE«XU) 

X(2,1)«X(2,LU>«U,0 
LL«U2 
DO  10  I ■2#LL 
J  ■  l-|*2 

A1«0.!»*CX(1,  I)*XI1,  J)| 
A2«Q,9*(X(2,l)»Xtl,J)) 
91*0»5*(*X(l«|)*X|l,j|) 
82*0»3*(*X(2#I>*X(2»J)) 

AI  ■  1-1 
Fl  >  F*ZI 
6l«C0SFOl) 

82**S INK (Fl ) 

SAVE-B1 

M1>B1^Q1-B2«G2 
B2»B2*Q1*SAVE*Q2 
X( 1, I  )bA1«B2 
X ( 2, I )■ A2+Bi 
X(1,J)  A1+B2 
10  X(2,J)«»A2+B1 

X(2,LL*1)«»X(2»LL*1) 

RETURN 

END 


oooonccooooooooooo 


SUBROUTINE  COQLH1.0R(N»X> 


THIS  COMPUTES  THE  HILBERT  TRANSFORM  OF  A  DATA  SERIES 
USING  THE  HYPER-RAPID  FOUKUR  TRANSFORM  ROUTINE  COOL 
THIS  program  ihanks  TO  jon  claerbout 


INPUTS  - 

N  *  LUG  (BASE  2)  OF  NUMBER  OF  DATA  POINTS 
REAL(X)  ■  DATA  SERIES  TO  BE  TRANSFORMED 
1  MAG ( X  )  •  0 

0UTPU1S  • 

REAL(X)  f  X  AGAIN 
I  MAGI  X )  •  HILBERT  TRANSFORM  OF  X 

THIS  CALLS  COUL 

DIMENSION  X ( 1 > 

TYPE  COMPLEX  X 
CALL  COOL (N*  X» "1 i 0 ) 

M  ■  2**N 
Ml  ■  M/2+2 
DO  1  1«M1»M 

i  mn  * 

X(l>  ■  ,s«xui 
X(Ml-l)  ■  ,S*X«M1.1» 

CALL  COOL<N»X(«1«0> 

RETURN 

END 


ouuuuuuouuuuouououuouuuu 


SUBROUTINE  C00UWQCN#X,SIQNf  A.tf) 


THIS  USES  COOL  TO  COMPUTE  THE  FOURIfcR  TRANSFORM  OT  TWO 
T I ME  SERIES  AT  ONCE 

INPUTS  - 

N  LOG  (BASb  2)  OF  NUMBER  OF  DATA.  POINTS 

X  A  COMPLEX  ARRAY  OF  DATA,  THE  FIRST  TIME  SERIES  IS  STORED 

IN  THE  REAL  PART  OF  X#  AND  THt  SECOND  IS  STORED  IN  THE 
IMAGINARY  PART  OF  X.  IN  OTHER  WORLS#  THE  TWO  SERIES  ARE  1 
MULTIPLEXED  IN  THE  ARRAY  X, 

SIGN  ■  -1«Q  FOR  DIRECT  TRANS?. iM,  THIS  SUBROUTINE 
HAS  NOT  BEEN  CHECKED-  OUT  f OH  TWO  INVERSE  TRANSFORMS 
AT  ONCE. 


OUTPUTS  • 

A  COMPLEX  FOURIER  TRANSFORM  OF  THE  FIRST  DATA  SERIES# 

I.E«#  THE  ONg  STORED  IN  THE  REAL  PART  OF  X 
B  FOURIER  TRANSFORM  OF  THE  SECOND  DATA  SERIES#  I.E.#  THE 
ONE  3T0RE0  IN  THE  IMAUINAHY  PART  OF  X, 

BOTH  TRANSFORMS  ARE  OF  LENGTH  2**CN-1»  ♦  1  ISEE  COOL  WRITEUPI 


DIMENSION  X(1)#A11),0(1) 

TYPE  COMPLEX  X#A#B,CONJG 
CALL  COOL(N#X«SlUN) 

ACT)  ■  ,5MXU|*CONJG<X<l»n 

BID  «  <U,,.,5>*<xtlj*CONJU(X(l)  )) 
M«2**N 

no  iu  k«2#h 

A(K)»0,t>*(X(K)  +  CUNJG(X(M  +  2-K)  I ) 

10  B(K)*(0«U*«0*9)*lX(K )»CUNJU(X(M*2*K * ) I 
RETURN 
END 


.  ■  t  rra  '  » 


*• 

'  2L  •<-■■■  - 


SUBRUUT  I Ng  CO0LVU'.vaX,X»U#M 
0 1  MENS  1  ON  FU)«X42«1) 

SINGLE-CHANNEL  CONVOLUTION  USING  COOL 


THIS  TAKES  FOURIgR  TRANSFORM  OF  UAT'i 
THfcM  I  OGETHfcH#  ANO  TRANSFORMS  BACK, 


INPUTS  . 


LX  LENUIH  OF  DATA 

LF  LENUIH  OF  FUTfcR 

►  MITER  COEFFICIENTS  DIMENSIONED  F < CF  >  IN  CALLING  PGM 

X  DATA#  DIMENSIONED  X f N>  IN  CALLING  PGM,  WHERE 

N  IS  THE  SMALLEST  NUMBER  WHICH  IS  A  POWER  OF  2  EXCEEDING 
<LF  *LX ) *2 

THE  SUBROUTINE  RETURNS  X  CONVOLVED  WITH  T,  OF  LENGTH 
lf+lx-i#  stored  CLOSE-PACKED  IN  X, 

23  SEPi’EMBEN  1966  UWMCC 


CHECK  LENGTH  RESTRICTION 

NX*LF*LX 
00  10  1*1,13 
N*2**l 

IMNX-NI  20,20,lu 
CONTINUE 

EHROR  RETURN  -  LENGTH  OF  FILTERED  RECJRD  WOULD  EXCEED  LIMIT 

LF«-LF 

RETURN 

NCOOL  *1 

ERASE  WORKING  SPA'-E  IN  X 

CALL  ERASE(2*;!"LX,XlLX«l>) 

MULTIPLEX  DATA  AND  FILTER  IN  X 

00  30  1*1, LX 
U*LX- I +1 
X(1,U)  «  X ( J > 

00  36  1*1, LX 
X 1 2, I )  •  O.Q 
DO  40  I*1#LF 
X ( 2, 1  )  *  Fill 

TRANSFORM  AND  FIDDLE 

FN*N 

CALL  COUL ( NCOOL*  X ,  <■! . 0 ) 

X(l, II  *  X ( 1# 1 ) *X ( 2, 1 1 /FN 


U  U  U  U  U  u 


$0 


60 


X(2,l)  ■  0,0 
N2«N/2 


00  50  !L"2,N2 

I  «5.  S  •  .y;  * -r 

1  M  #  *FN J  *  "  IU+Z,**2*X<i»n.)**2*X<2,lL)*#2)/ 

x  c  1 , 1 L » - 1 


X < 2# N2*l ) «0 ( Q 
N22*N2*2 


00  60  IL*N22#N 
Ml,  !LMX<l#N«IL*a$ 
x < 2,  IL >•  ~X<2,N» JL+2) 


transform  back 
CALL  COOL <NCOOL#*, +1,0) 

OLOSB-PACK  FILTERED  DATA  in  X 

00  70  IH1»N X 
*111  ■  XU, I) 


return 

6NU 


ooo 
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SUBROUTINE  FT2DCU0L<X#N#M#SIUN| ) 
dimension  *<n#mi 
TYPE  COMPLEX  X 

2  DIMENSION  FUURlER  TRANSf OHM  USING  COOL 

NC00l*L0UF< FLOM*  |N)  >/LOUF  ( 2 . 0  )  +X ,  Ufc-4 
MCOOL»LOUF|FLOAT* f  M  > )/LOQF<2tOUl,Ofc-6 
SCALfcN«l*0/SQRTMFLOATFCNI ) 

9CALfcM*l , 0/SQRTF (RLOaTFIM) ) 

00  1  IM«1,M 

CALL  COOLCNCOOL* X*l# IM) jSJUNI) 

00  1  IN*1»N 

1  X( IN# IM)>X( IN, |MI*$CALEN 
CALL  MATHA69(X#N#M#X) 

00  2  IN*1#N 

index«i«(in*d*m 

CALL  COOL<MOOOL#X| INDEX# 1)#S1QNI) 

CALL  SCALE! SC ALEMiH#X( IN0EX#1) ) 

2  CONTINUE 

CALL  MATKA63<X»M#N#X) 

RETURN 

KNU 


non 


subroutine  ft3ooooux#n#n,l,si«nI) 

0 1  MENS  I  ON  X(N«M#L) 

TYPE  COMPLEX  X 

3  DIMENSION  FOURIER  TRANSFORM  USING  COOL 

LCOOL*LOOF(FLOA7F  ILM/LQQFI2.0U1. OE-6 
SCALtL"!*  O/SORTF <FLOATF(L ) > 

OO  1  IL*1«L 

CALL  FT2UC00LIXUtl#!L)«NfM(SIGNll 

1  CONTINUE 

CALL  MATKA63<X*N*M,L,X) 

DO  2  IN"1,N 
OO  2  IM"1#M 

lNDEXal+( lN«l)*L^I 1M«11*L*N 

CALL  COOL(LCOOL*X( INDEX* 1*1) *  SIGN J ) 

CALL  SCALE! 8CALEL|L*X(IN0EX« 1*1)) 

2  CONTINUE 

CALL  MATRA63<X*L«N*M«X! 

RETURN 

END 


SUdHUUTINfc  MAJHA63( A#N,M*8» 

U 1 MfcNS I  UN  A<2,1)*B(2,1> 

C 

C  MATRIX  THANbMUSE  ON  COMMlbX  AKRAYS 

C 

MASK  1*0 OOOOOOQOOUQOOOlB 
MASK2«7;/777/77//77776B 
NM*N*M 

00  10  i*i,nm 

d(l, I )*A(1, I ) ,OR,MASKl 
10  B(2« I )*A<2, I > 

Jr  *0 

ASSIGN  3U  TO  KSHH 
00  100  I *1#  NM 
GO  TU  KSwH, <3Q,!>0| 

30  JF«JM1 

LL*B(1#JF),AND*MASK1 
I F  < ll ) 30 *  30 , AQ 
40  JO*  JF  <*1 

ASSIGN  SO  TO  KSMH 
?6MPdl*B ( 1, JF  ) 

TEMPB2*BC2, JF  ) 

50  Jl*JO/N*XMOl)F| JU«N)*M*1 
TEMPA1=B(1,J1) 

TEMPA2*B(2, Jl) 

B<1*J1)*TEMPB1.AND,MASK2 

d<2, J1)*TEMP82 

I EMPdl*  TfcMPAl 

rEMPB2rTfcMPA2 

J0«J1-1 

IF  CJ1-JF >60,60,100 
60  ASSIGN  30  TO  KSwM 
100  CONTINUE 
HfclUHN 
fcNU 
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SUBROUTINE  SPfcCINUMt IT,UT,KT»X,L#LF,S> 
DIMENSION  X<2, 11, 9(2, 1 » 

DIMENSION  X(2,N«N,LF),X(2»;?*LX*2),S(2,LF) 


SPECTRAL  MATHIX  FOR  TAPED  data 

DATA  MUST  BE  A  POWER  OF  TWO  |N  LENGTH  AND  ON  TAPE  IT 
FORMA!.  SPECTRAL  MATRIX  IS  RETURNED  AS  A  F6J  COMPLEX 
IN  X < 

1 1  *  input  subset  tape 
jt-schatch 

K  T  ■•SCRATCH 

X-WORKING  AHHAy  and  returned  spectral  matrix 
L-skUMBER  of  times  to  smooth 
lf-returneo  length  of  spectral  estimates 

S-WOHKING  AHHAy 

PROGRAM  TOO  COMPLICATED  to  DESCRIBE,.. 

rewind  ir 
REMIND  OF 
REWIND  KT 
READUT)LOST,N»LX 
LX2-2-LX 

NCOOL-LOGFc FLOATS  ILX) J/LOGF  l*,U)n*PE*6 

NSCI-N-N 

LF«LX/2«*L4l 

LX2P2*LX242 

LX4-A-LX 

LX2P2T2»2*LX2P2 

LXP1-LX*! 

LX2P<J*LX2+3 

IDC-0 

LF2-2-LF 

WRITE(JT)L0ST,N,LX2P2 

WRITE(KT)LOST,N.LX2P? 

DO  10  IN-1, N 

CALL  ERASE (LX4»  X I 

READCITMX(l,M»,Mwl,LX) 

Call  COOL (NCOOL  +  lfx,<*l,Q) 

WRITE  (JTMX(M),M-1#LX2P2) 

WRITE(KTMX(M|,M-X,LX2P2» 

10  CONTINUE 
END  ► ILE  JT 
END  F ILE  KT 

rewind  ir 
rewind  jt 

REWIND  ki 
00  1  IN-1, N 
I ND- I N*1 

CALL  SK  l  PREC  < I N»  K T ! 

READ(KTHX<M),M*laLX?P2) 

CALL  D0TEMIX»X#LAP1,XCLX2P4> » 


IN  SUBSET 
MATRIX 


CALL  SM0UTH<X<LX*p3),LXPl»L) 

CALL  DISC63< IDC*1,X(1  X2PJ)»LF2> 

IOC* I DC*V 
UO  6U  JN“ IND*  N 

HEA0<KTMX<M>,M*Ln2P3*LX2P*T*> 

CALL  OOTEMf  X*  X(LX2P3)*LXP1*  X<LX2P3>  > 
CALL  SM00TH(X(LX2P3I*LXP1»L) 

CALL  0ISC63(IDC*1,X<LX2PJ)#L*’2> 

IDC* 1 DC*V 
60  CONTINUE 
REWIND  KT 
ISAVfc*KT 
KT*JT 
JT" I  SAVE 

1  continue 

IDC*U 

00  25  IN*i,N 
I ND* I N*1 

CALL  0ISC63(IDC*UiS»LF2) 

I  DC* I  DC* V 
INOfcX«IN*( IN-1»*N 
DO  26  IL«1,LF 
X(l,  INDEX)«S<i»  ID 
X(2,INDEX)»S<2»1L) 

INDEX»INUEX*NSU 
26  CONTINUE 

00  it  JN*lND.N 

CALL  0ISC63(10C*U,S*rF2) 

IDC»IUC*V 

IN06X1«1N*( JN,1I*N 
INDEX2*JN«(1N>.1I*N 
00  2b  IL»1,LF 
X(l* INDEX1)*S(1» 1L) 

X  ( 2  #  INDEX  1)*S(2*  ID 
X(l.  I  NDEX2 ) *S ( 1* 1L) 
X(2»IUDEX2)**8(2*IL) 

iNDEXl^INDEXltNbU 

INDEX2«i|NDEX2^NSu 

26  CONTINUE 

27  CONTINUE 
25  CONTINUE 

RETURN 

END 


non 


subroutine  smquim|x#i.en<jth#l> 


THl  i  HANNING  NQUTINE  THANKS  TO  U  CLAERBOUT 


DIMENSION  KC2#LfeNQTH} 

Lf  *LfcNGTH 
I.FM1-LF-1 
DO  1  IL*1»L 

<  < 1* 1 )  =  U.5*X(1»1)*0.5*X(1*2) 
X|2,l)«o»0 

X  (1#LF  '*0,5#X(1#LF)*o,5*XU,LF-1) 

X ( 2#LF ) »0 , 0 

IND*2 


DO  2  JLB  i#LFMl*2 

X(2*JL)«U,25*X(2*  JL»1)*0.5*X<2#  Jl >*0 ,25*X ( 2. JL+l > 

*n.'w;.;?n3u'JL'U*0'S',l'1'J1-l*0,25*xu'JL*i' 

*(2» IND)»X(2# JL) 

ind«ind*i 

CONTINUE 


X  C 1#  I  NO)  *X  ( 1#  l,F  | 
X ( 2# 1  NO ) aX ( 2#  LF  ) 
IF  *LF /2*1 

lfmi-lf-i 

1  CONTINUE 
RETURN 
End 


l. 


SUBROUTINE  DOTbM < X, Y, L, 1 1 
DIMENSION  x C 2#L) * V 1 2*L > * 2 (2* L ) 
DO  1  IL*1»  L 


SAVER*XU#  IL)AY(1*IL)*X(2«1L>*Y(2,IL) 
SAVE  I  «X  (1*  Il>*Y(2#  IU-XI2#  JLT*Y(J,IL) 
2 ( i#  I L )  "SAVER 
2(2. 1 1  > "SAVE  1 
ML  TURN 


END 


onoooooooooooon 


SUBROUTINE  Ul  SOO-if  I  BLOCK#  IBB  I  TCH^X.N) 

DIMENSION  y ( N  > 

THIS  IS  THg  BUt  DISC  DRIVER  ROUTINE  WRITTEN  JN  CODAP-1 
IT  TRANsFfcRS  Words  BETWEEN  CORE  AND  THE  DISC 
I  BLOCK  is  the  DISC  BLOCK  (32  WORDS)  ADDRESS 
1SWITCH  CONTHUtS  READING  ANU  WRITING 

I SW I T QH» U  GIVES  A  HEAD  FROM  THE  DISC 
I SW I T QH»1  GIVES  A  WRITE  ON  THE  DISC 
X  IS  THE  CORE  ADDRESS 
N  IS  THE  NUHdfcR  OF  WORDS  10  TRANSFER 


THIS  ROUTlNb  «UST  BE  SUPPLIED  BY  THE  USER  OR  INCLUDED  IN  BINARY 


RETUHN 

END 


SUBROUTINE  ERABE(N.X) 
DIMENSION  x(N) 

C 

C  ERASE  N  WORDS  IN  X 

C 

00  1  I »1»  N 
1  X(I)*0,0 
RETURN 
ENO 


SUBROUTINE  SK{PHtC(N, ITAPS) 

SKIP  N  LOGICAL  RECORDS  UN  TApfc  { TAPE 

DO  1  I  *  1 » N 
1  READ ( IT  APE ) LOST 
RETURN 
END 


APPENDIX  C  -  PROGRAM  WRITE-UPS 
FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA , VIRGINIA 

digita:  computing  section 

A.  IDENTIFICATION 

Title;  Hyper-Rapid  Specialized  Cooley-Tukey  Fourier  Trans¬ 
form  (direct  only) 

COOP  Identification;  G612-COOL 

Category:  Fourier  Transform 

Programers:  J.  F.  Claerbout,  D.  W.  McCowan,  J.  L.  Gibson, 
and  E.  A.  Flinn 

Date;  26  February  1966 

B.  PURPOSE 

To  compute  the  Fourier  series  expansion  of  a  real-or 
complex-valued  data  series,  or  the  data  series  from  the 
complex-valued  Fourier  series  expansion. 

C.  USAGE 

1 .  Operational  Procedure  and  Parameters: 

This  is  a  CODAP  subroutine  with  a  FORTRAN-63 
calling  sequence  CALL  COOL  (N,  X,  SIGN)  .  X  is  a  com¬ 
plex  array  used  for  the  data  series  and  the  transform; 

N 

the  number  of  elements  of  X  is  L  =  2  ;  SIGN  =  -1.0  for 
a  direct  Fourier  transform,  and  +1.0  for  an  inverse 
Fourier  transform  (but  see  below  for  arrangement  of 
data) . 

* 

For  the  direct  transform:  on  input  the  real  part 
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Of  X  contains  the  data  serins  and  the  imaginary  part  of  x 

”  Zer°-  °n  retUrn'  the  F°“ier  cosine  series  expansion 
is  in  the  real  part  of  X,  and  the  Fourier  sine  series 

expansion  is  in  the  imaginary  part  of  X.  Each  contians 

only  2  +  i  nonredundant  points;  the  cosine  expansion 

is  symmetric  about  point  number  2N_1  +  i  and  the  sine 

transform  is  antisymmetric  about  this  point. 

For  example;  N  =  3  and  data  =  (0.,  1.,  0.,  o.,  0., 

°-'  #-  Re(X>  "  (0-'  0.,  0.,  0.,  0.,  0.); 

In(X)  ‘  (0"  °-'  *•'  "•>  «•.  0.,  0.)  on  imput.  on 

return,  Re(x>  -  (l.ooo.  .7071,  0.,  -.7071,  -i.ooo,  -.7071, 

°-'  -7071,;  Im<X)  =  <°->  --7071,  -l.ooo,  -.7071,  0., 

•7071,  1.000,  .7071).  Point  number  1  corresponds  to  xero 
frequency;  point  number  5  corresponds  to  it. 

.  POr  inverse  transform:  the  cosine  and  sine  series 
must  be  folded  over  about  point  number  21,':1  +  1  before 
calling  COOL  with  SIGN  =  +l.o. 

There  is  a  scale  factor  of  2_N  which  COOL  does  not 

apply.  The  user  can  choose  to  apply  the  scale  factor 

Githsr  to  th©  dirsct  o r  f ^  _  • 

®Ct  or  to  the  inverse  transform,  or  to 

apply  a  factor  of  2^  to  both.  For  example,  if  cool 

were  called  with  the  transform  example  above,  the  result 

would  be  Re (X,  -  (o.,  8.,  0.,  0.,  0.,  0.,  0.)  and 

Im(x)  =  (0-'  °-'  °w  o.,  0.,  0.,  0.,  o.). 
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2.  Space  Required:  Approximately  20010  exclusive  of  X. 

The  largest  series  that  can  be  transformed  in  a  32K 
core  machine  is  8K. 

3.  Temporary  Storage  Required;  None.  Other  versions  of 
this  program  have  an  auxiliary  storage  for  the  cosine 
table  and/or  a  table  of  bit-reversed  numbers.  COOL 
computes  its  sines  and  cosines  as  it  goes,  and  uses  an 
algorithm  due  to  J.  F.  Claerbout  for  calculating  the 
bit-reversed  numbers. 

4.  Frintout ;  None. 

5.  Error  Printouts:  None, 

6-  Error  Stops:  None. 

7.  Input  ard  Output  Tape  Mountings:  Not  Applicable 

8.  Input  and  Output  Formats;  Not  Applicable. 

9.  Selective  Jumps  and  Stops;  None. 

N 

10.  Timing:  Time  is  proportional  to  N'2  .  Transforming 
8192  on  the  CDC  I6't4-B  requires  25.0  seconds. 

11.  Accuracy:  Calling  COOL  returns  the  original  to  about 
nine  decimal  places. 

12.  Cautions  to  User:  See  Operational  Procedure  above. 

13.  Configuration:  Standard  COOP. 

14.  References:  J.  W.  Cooley,  1964  "Harm  -  Harmonic  Analy¬ 
sis;  Calculation  of  Complex  Fourier  Series":  IBM  Watson 
Research  Center  Yorktown  Height,  New  York. 
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J.  W.  Cocley  and  J.  W.  Tukey,  1965,  An  Algorithm  for 
the  Machine  Calculation  of  Complex  Fourier  Series: 
Math,  of  Comp.,  Vol.  19,  pp.  297-301. 


Writeups  of  the  following  SDL  programs: 
COOLTWO:  Does  two  Fourier  transforms  at  once. 

FT3DCOOL:  Three-dimensional  Fourier  transform 

D.  METHOD 

N 

Given  a  time  series  X(I),  1,  L  (where  L  =  2  )  assumed 
to  be  periodic  outside  the  given  range,  COOL  constructs 

N-l 

Y (K)  =  SUM  X(J)*WJK  K  -  0,  L  -  1 

J=0 

where  W  ®,  exp  (-2rri/L)  for  time-frequency  transform,  and 
W  *=  exp  (+2rri/L)  for  frequency-time  transform.  The  algo- 

,  ,  #  JJ 

nthm  is  efficient,  requiring  N’2  multiplications  rather 
2N 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

IDENTIFICATION 

Title;  Multichannel  convolution  in  the  frequency  domain, 
for  taped  data., 

COOP  Identifications  UES  G620  COOLCON 
Category:  G6  Time  Series  Analysis 

Programer :  D.  W.  McCowan 

Date;  22  September  1966 
PURPOSE 

This  subroutine  convolves  data  channels  on  the  input 
subset  tape  with  a  multichannel  filter  stored  in  core, 
working  entirely  in  the  frequency  domain.  The  result  is 
written  in  subset  format  on  another  tape. 

USAGE 

1.  Operational  Procedures  This  is  a  FORTRAN-63  subroutine 
with  calling  sequence; 

CALL  COOLCON  (INT ,  IOT,  L,  F,  X). 

2 .  Parameters; 

INT  is  the  number  of  the  input  tape  unit. 

IOT  is  the  number  of  the  output  tape  unit. 

L  is  the  number  of  points  in  the  filter  (see  restriction 
below)  . 
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F  is  the  multichannel  filter,  dimensioned  F(N,L)  in  the 
calling  program,  Where  N  is  the  number  of  channels  on  the 
input  subset  tape. 

X  is  a  working  array,  dimensioned  X(2,IT)  is  the  calling 
program,  where  IT  is  the  l(east  power  of  2  such  that 
IT 

2  >  L  +  LX 

where  LX  is  the  number  of  data  points  in  the  input 
channels . 

Restriction  on  length  of  data  and  length  of  filter: 

13 

LX  +  L  must  not  be  greater  than  2  (8K) . 

3.  Space  Required;  Very  little  in  addition  to  arrays. 

IT 

4.  Temporary  Storage  Required:  2*2  working  space,  plus 
127^  for  the  subset  tape  label. 

5 •  Printout  s  None . 

13 

6.  Error  Printouts:  If  L+LX>2  ,  these  numbers  are  printed 
with  an  error  message. 

13 

7.  Error  Stops:  If  L+LX>2  ,  the  subroutine  stops  the 
calling  program. 

8 .  Input  and  Output  Tape  Mountings:  See  Parameters  above . 

9.  Input  and  Output  Formats:  Compatible  with  UES  Subset 
(See  Writeup)  . 

10.  Selective  Jump  and  Stop  Settings;  None. 

11.  Timing:  Dominated  by  two  Fourier  transforms  using  COOL 
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The  length  of  trans- 


for  each  channel  to  be  filtered. 

XT 

form  is  2  (See  Writeup  o£  COOL) . 

Accuracy:  This  yields  the  same  numbers,  to  ten  decimal 
places,  which  would  be  computed  by  convolving  the 
filter  and  data  series  in  the  usual  way. 

Cautions  to  User:  None . 

Configuration:  Standard  COOP. 

References:  Writeups  of  UES  Gfal2  COOL,  UES  Z24  SUBSET, 
and  UES  G617  COOLER. 

METHOD 

Fox  each  channel  to  be  filtered,  the  subroutine  erases 
-IT+1  1 

locations  of  X,  and  multiplexes  the  filter  and  the 
data  channel  in  X,  starting  at  the  beginning.  Note  that  as 
far  as  COOL  is  concerned,  X  is  a  complex  array  with  data  in 
the  real  part  and  filter  in  the  imaginary  part.  COOL  is 
called,  and  the  logic  of  COOLER  (q.v.)  i»  used  to  form  the 
Fourier  transform  of  the  filtered  channel  in  X.  COOL  is 
called  again  to  get  back  to  the  time  domain,  and  the  filter¬ 
ed  channel  is  written  on  the  output  tape. 

The  subset  label  is  copied  from  the  input  tape  to  the 
output  tape  at  the  beginning  of  the  subroutine. 


12. 

13. 

14. 

15. 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

IDENTIFICATION 

Tijble:  Hyper-Rapid  Specialized  Cooley-Tukey  Fourier  Trans¬ 

form  (direct  only) 

COOP  Identification:  G617-COOLER 

Category:  Fourier  Transform 

Programer:  J.  F.  Claerbout 

Date:  27  July  1966 

PURPOSE 

To  compute  the  Fourier  series  expansion  of  a  real-valued 
time  series. 

USAGE 

1*  Operational  Procedure;  This  is  a  FORTRAN-63  subroutine, 
with  calling  sequence  CALL  COOLER(N,X).  This  subroutine 
calls  COOL. 

2*  Parameters :  On  input,  X  is  a  real-valued  time  series 
containing  LX  points,  where  LX  =  2N,  N  is  restricted 
to  be  14  or  less.  On  return,  X  contains  *sLX+l  complex 
points  of  the  Fourier  transform  of  the  data,  with  the 
real  and  imaginary  parts  multiplexed  together  -  i„  e., 
on  return  X  can  be  thought  of  as  a  complex  array,  with 
the  cosine  transform  in  the  real  part  and  the  sine 
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transform  in  the  imaginary  part. 

X  must  be  dimensioned  at  least  LX+2  in  the 
calling  program.  (i.e.f  JsLX+l  complex  points) 

3.  Space  Required:  Very  little. 

4.  Temporary  Storage  Required:  None. 

5 .  Printout:  None . 

6.  Error  Printouts:  None. 

7 .  Error  Stops:  None . 

8 .  Input  and  Output  Tape  Mountings:  Not  Applicable . 

9 .  Input  and  Output  Formats:  None . 

10.  Selective  Jumps  and  Stops:  None. 

N 

11.  Timing:  Time  is  proportional  to  N.2  ;  transforming 
16384  points  on  the  CDC  1604-B  requires  45.0  seconds. 

12.  Accuracy:  About  nine  decimal  places. 

13.  Cautions  to  User:  On  return,  the  real  and  imaginary 

parts  of  the  transform  are  multiplexed  together.  X 
must  be  dimensioned  at  least  LX+2  in  the  calling  pro¬ 
gram,  not  LX.  This  subroutine  will  not  do  an  inverse 
transform. 

14.  References:  Writeup  of  UES  G612  COOL. 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

IDENTIFICATION 

Title:  Hilbert  transform  of  periodic  data 
COOP  Identification:  UES  G619  COOLHLBR 
Category:  G6  Time  Series  Analysis 
Programer:  E.  A.  Flinn  and  J.  F.  Claerbout 

Date:  23  September  1966 
PURPOSE 

To  compute  the  Hilbert  transform  (quadrature  function) 
a  i-ime  series.  Since  COOL  is  used,  the  time  series  is 
assumed  to  be  periodic  outside  the  range  of  definition. 

USUAGE 

1*  Operational  Procedure:  This  is  a  FORTRAN-63  subroutine, 
with  calling  sequence:  CALL  COOLHLBR(N,X) .  This  subroutine 
calls  COOL. 

2*  Parameters :  N  is  the  log  (base  2)  of  the  number  of 

data  points.  X  is  the  data,  dimensioned  at  least  2N  in  the 

calling  program,  and  type  complex  there. 

On  input,  the  real  data  series  must  be  stored  in  the 

real  part  of  X,  and  the  imaginary  part  must  be  zero. 

On  return,  the  real  data  series  is  stored  in  the  real 

1 

part  of  scaled  up  by  2  "  .  The  Hilbert  transform  is  stored  in 
the  imaginary  part  of  X,  also  scaled  up  by  2N_1. 
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Space  Required;  Very  little  in  addition  to  the  array 

N+l 

for  data,  which  requires  2  locations  in  the  calling 
program . 

Temporary  Storage  Required;  None 
Printout :  None 
Error  Printouts;  None 
Error  Stops;  None 

Input  and  Output  Tape  Mountings:  Not  Applicable 

Input  and  Output  Format ss  Not  Applicable 

Selective  Jumps  and  Stops;  None 

Timing:  Dominated  by  two  calls  to  COOL 

Accuracy:  The  data  is  returned  correct  to  ten  decimal 

places . 

Cautions  to  User:  The  data  must  bn  arranged  as  under 
( 2 )  above . 

Notice  that  as  far  as  this  subroutine  is  concerned 
the  data  is  periodic  outside  the  range  of  definition .  End 
effects  rftay  cause  answers  which  the  user  does  not  expect 
For  example,  if  the  input  is  a  pure  sine  wave,  the  user 
expects  the  quadrature  to  be  a  pure  cosine.  Using  this 
subroutine,  this  turns  out  to  be  the  case  only  if  the  data 
series  contains  an  integral  number  of  cycles. 

References :  Writeup  of  UES  G612  COOL. 


D .  METHOD 


The  Hilbert  transfer™  of  a  function  has  a  Fourier  trans¬ 
form  which  is  (-1}*5  times  -the  Fourier -transform  of  the  function. 
COOL  returns  the  real  and  imaginary  parts  of  the  Fourier  transform 
of  a  function  calculated  from  zero  to  2n,  so  that  the  real  part 

ia  symmetric  about  the  middle  and  the  imaginary  part  is  anti- 
symmetric. 

If  the  Fourier  transform  of  the  function  is  A+iB,  the  Fourier 
transform  of  the  Hilbert  transform  is  -B+iA.  All  COOLHLBR  does 
is  erase  the  second  half  of  the  Fourier  transform  (the  part 

from  tt  to  .tt)  ,  half-weight  the  end  points,  and  call  CCJL  again 
to  transform  back  to  the  time  domain. 

The  scale  factor  2N-1  comes  from  the  fact  that  COOL  gives  J 

the  unnormalized  transform.  I 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A.  IDENTIFICATION 

Title:  Fourier  Transform  of  Two  Data  Series  Simultaneously 
COOP  Identification:  COOLTWO 
Category:  G6  Time  Series  Analysis 

f 

Programer:  E.  A.  Flinn 

Date:  10  June  1966 

B.  PURPOSE 

To  compute  the  Fourier  series  expansion,  using  COOL 
(q.v.),  of  two  data  series  simultaneously. 

C.  USAGE 

1*  Operational  Procedure:  This  is  a  FORTRAN-63  subroutine 
with  calling  sequence. 

CALL  COOLTWO  (N,  X,  SIGN,  A,  B)  . 

2 .  Parameters: 

N  is  the  log  (base  2)  of  the  number  of  elements  in  X; 

X  contains  the  two  data  series,  multiplexed  in  one  complex 
array,  so  that  Re(X)  contains  one  series  and  Im(X)  contains 
the  other . 

SIGN  =  -1.0  .  The  program  has  not  yet  been  checked 
out  for  inverse  transformation; 

A  is  the  complex  (cosine  and  sine)  transform  of  the  data 
series  stored  in  the  real  part  of  X; 


-  C-l  3 


B  is  the  complex  Fourier  transform  of  4-he  data  series  stored 

■i 

in  the  imaginary  part  of  X; 

A  and  B  are  both  of  length  2** (N  -  1}  +  1. 

3.  Space  Required  a  about  70.^  excluding  arrays  , 

4.  Temporary  Storage  Requirements s  None 

5 .  Printouts;  None 

6.  Error  Printouts;  None 

7 .  Error  Stops:  None 

8 .  Input  and  Output  Tape  Mountings;  None 

9.  Input  and  Output  Formats;  Not  Applicable 

10.  Selective  Jump  and  Stop  Settings;  Not  Applicable 

N 

11.  Timings  Timing  is  proportional  to  N-2  ;  transforming 
8192  data  points  on  the  CDC  1604-B  requires  25.0  seconds. 

12.  Accuracy;  Same  as  COOL. 

13.  Cautions  to  Users  This  program  has  not  been  checked  out 

for  inverse  transformation.  This  program  does  not 

-N 

apply  the  scale  factor  0  t  since  some  users  may  wish 
to  apply  the  scale  factor  to  the  inverse,  rather  than 
the  direct  transform.  The  number  of  data  points  must 
be  a  power  of  2 . 

14.  Configuration s  Standard  COOP 

15.  References Writeup  of  UES  G612  COOL 
D.  METHOD 

The  method  is  due  to  J.  W.  Cooley  (see  Reference  2  in 
main  body  of  this  report.) 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA ,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A,  identification 

Title:  Fast  convolution  of  two  time  series  using  CQQE'o 

COOP  Identification:  UES  COOLVOLV 

Category:  Time  Series  Analysis 

Programer s  E.  A.  Flinn  and  D .  W.  McGowan 

Date:  23  September  1966 

B.  PURPOSE 

/ 

To  form  the  convolution  of  two  time  series,  not  by  the 
usual  polynomial  multiplication  algorithm,  but  by  forming  the 
two  Fourier  transforms  (using  COOL),  multiplying  them  together, 
and  transforming  back  to  the  time  domain.  This  is  faster  than 
the  usual  procedure  when 

LX-LF  »  4  (2N  +  1)  (LX  +  LF ) 

where  LX  is  the  data  series  length,  LF  is  the  filter  impulse 
response  length,  and  N  is  the  log  (base  2)  of  LX  +  LF. 

C„  USUAGE 

1.  Operational  Procedures  This  is  a  FORTRAN-63  subroutine, 
with  calling  sequences 

CALL  COOLVOLV ( LX , X , LF , F ) 

2 .  Parameters : 

X  is  the  data  series  to  be  convolved,  dimensioned  at  least 
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2  in  the  calling  program,  where  2°  is  the  smallest 
power  of  two  larger  than  LX  +  LF« 

LX  is  the  length  of  the  data  series  to  be  convolved . 

F  is  the  filter  to  be  convolved  with  X. 

LF  is  the  length  of  the  filter. 

3.  Space  Required;  30°io  plus  arraYs • 

4.  Temporary  Locations  Required;  None  beyond  filling  out  X 
to  the  first  power  of  two  greater  the  LX  +  LF. 

5.  Alarms  or  Special  Printout;  None 

13 

6.  Error  Returns;  If  LX  +  LF  >  2  ,  LF  is  replaced  by  -LF 

and  control  is  returned  to  the  calling  program. 

7 .  Error  Stops;  None 

8.  Tape  Mountings;  None 

9.  Formats:  None  a 

10.  Jump  and  Stop  Settings:  None 

11.  Timings  Dominated  by  two  calls  to  COOL  for  LX  +  LF 
points  each  time. 

12.  Accuracy:  Gives  the  same  results  as  polynomial  multi¬ 
plication  to  ten  decimal  places. 

13.  Cautions:  None 

14.  Configuration:  Standard  COOP 

15.  References:  Writeups  of  COOL,  C00LC0N,  AND  COOLER 
METHOD 


The  same  method  is  used  as  used  in  C00LC0N. 


SEISMIC  DATA  LABORATORY 
ALEXANDRIA, VIRGINIA 

DIGITAL  COMPUTING  SECTION 

IDENTIFICATION 

Title;  Two  and  Three  Dimensional  Fourier  Transform  Package 

COOP  Identifications  G615  FT2DCOOL ,  FT3DC00L 

Category;  G6  Time  Series  Analysis 

Proqramer ;  D.  W.  McCowan 

Date :  20  April  1966 

PURPOSE 

The  subroutines  in  this  package  compute  two  and  three 
dimensional  Fourier  transforms.  Their  names  ares  FT2DC00L , 
FT3DC00L ,  COOL,  MATRA63 ,  and  SCALE.  As  with  COOL,  the 
dimensions  on  the  data  must  be  a  power  of  two. 

USAGE 

1 •  Calling  Sequences 

CALL  FT2DCOOL  (X,N,M,  SIGNI) 

and 

CALL  FT3DCOOL  (X,N,M,L,  SIGNI) 

2 .  Arguments  s 

X,  the  complex  array  in  which  the  data  is  supplied  and 
in  which  the  Fourier  transform  is  returned.  If  real 
data  is  supplied,  it  must  be  put  into  the  real  part  of 
X  and  the  imaginary  part  must  be  erased. 


N,M,L,  the  dimensions  of  X„  Each  of  these  numbers  must  be 
a  power  of  two.  The  number  of  complex  points  in  che  Fourier 
transform  will  be  N/2  +  1,  M/2  +  1,  and  L/2  +1  in  each  direct¬ 
ion  . 

SIGNI ,  a  switch  determining  the  type  of  transform  to  be  per¬ 
formed.  SIGNI  =  -1,0  gives  a  direct  transform  (time  to  fre¬ 
quency),  and  SIGNI  =  +1.0  gives  the  inverse. 

3.  Space  Required;  500  locations, 

4 .  Temporary  Storages  None 

5 •  Alarms  and  Printouts;  None 

6 .  Error  Returns s  None 

7 .  Error  Stops;  None 

8 .  Tape  Mountings;  None 

9 .  Formats;  None 

10.  Jumps  arid  Stop  Settings;  None 

11.  Time  Required;  Three-dimensional  Fourier  transforms  require 
NM  +  NL  +  ML  one-dimensional  Fourier  transforms.  Two-dimen¬ 
sional  Fourier  transforms  require  N  +  M  one-dimensional 
Fourier  transforms.  For  the  timing  of  one-dimensional  Fourier 
transforms,  see  References. 

Accuracy;  Same  as  COOL 
Cautions  to  Users;  None 
Equipment  Configurations  Standard  COOP 
References;  Writeup  of  UES  G612  COOL  3/30/66 


12. 

13. 

14. 

15. 
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Do  METHOD 

The  direct  2  and  3-dimensional  Fourier  transforms  are 
defined  ass 
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Where  =  exp(2rri/N);  W2  =  exp(2ni/N);  W^  =  exp(2ni/L) 

The  two-dimensional  transform  is  broken  up  ipto  N  +  M 
one-dimensional  transforms  and  the  three-dimensional  transform 
is  broken  up  into  L  two-dimensional  transforms  and  NM  one-dimen¬ 
sional  transforms.. 
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SEISMIC  DATA  LABORATORY 
ALEXAND  RI A ,  V I RG  IN  I A 

DIGITAL  COMPUTING  SECTION 

IDENTIFICATION 

Titles  Spectral  Matrix  Estimates 
COOP  Identifications  G618  SPECTRUM 
Categorys  Time  Series  Analysis 
Programer s  D.  W.  McCowan 
Dates  10  July  1966 

m 

.PURPOSE 

This  is  a  package  of  three  FORTRAN-63  subroutines  for 
computing  an  estimate  of  the  spectral  matrix  for  N  channels 
of  data  stored  on  magnetic  tape.  It  uses  the  hyper  “-rapid 
Fourier  transform  routine  COOL,  and  makes  use  of  two  tapes 
and  the  disc  to  cut  running  time  to  a  minimum.  The  names 
of  the  three  routines  in  the  package  ares  SPECTRUM,  DOTEM, 
and  SMOOTH.  In  addition  to  these,  three  more  subroutines 
are  assumed  to  be  on  the  system  tape;  they  ares  COOL, 
SKIPREC,  and  ERASE .  Since  all  other  routines  are  called 
internally  by  SPECTRUM,  only  the  calling  sequence  for  it 
will  be  given. 

USAGE 

1 »  Calling  Sequences 

Call  SPECTRUM.  (IT,  JT,  KT,  S,  NS,  LF ,  X) 
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2 .  Arguments; 


IT,  the  input  subset  tape  number  on  which  the  N  channels 
of  data  are  written .  The  length  of  each  channel  must 
be  exactly  a  power  of  two. 

JT ,  the  number  of  a  scratch  tape „ 

KT,  the  number  of  a  scratch  tape, 

S,  a  triply  subscripted  FORTRAN- 63  complex  array  used  both 
for  internal  manipulation  and  to  return  the  computed 
spectral  matrix  as  a  N  by  N  by  LF  complex  array  with 
subscripts  varying  in  that  order „  Here  N  is  the  number 
of  channels  read  from  the  input  tape  label  and  LF  is  the 
smoothed  length  of  each  spectral  estimate .  This  array 
must  also  be  4*LX+4  locations  in  length,  since  it  is 
also  used  for  internal  computations „  LX  is  the  length 
of  the  input  data  channels  read  from  the  input  tape 
label.  Remembering  that  there  are  two  locations  used 
for  each  complex  number,  the  total  dimensions  on  S  in 
the  main  program  must  be  2*N*N*LF  or  4*LX+4,  whichever 
is  the  rger •  i-s  usually  convenient  to  dimension  it 

as  complex  N  by  N  by  L  F6,3  array  in  order  to  facilitate 
use  -  L  here  is  a  number  chosen  so  that  S  will  be  large 
enough  as  described  above. 

NS,  the  number  of  times  to  apply  the.  hanning  smoothing 
operation  to  the  original  estimates. 
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3. 

4. 

5. 

6. 

7. 

8. 
9. 

10. 

11. 

12. 

13. 


14. 

15. 


LP ,  the  returned  length  of  the  spectral  estimates .  This  is 
computed  from  the  formula: 

LF  =  (LX/(2**NS)  +  1 
LF  must  not  be  larger  than  129. 

X,  an  array  used  for  internal  manipulation,  containing  at 
least  2*LF  locations. 

Space  Required:  502  locations 
Temporary  Locations:  None 
Alarms  or  Special  Printout:  None 
Error  Returns:  None 

Error  s tops:  The  subroutines  stop  if  length  of  data  series 
13 

exceeds  2  . 

Tape  Mountings:  See  Arguments 

Input  and  Output  Formats:  See  Arguments 

Jump  Settings:  None 

Time  Required:  A  10-channel,  4096-point,  NS  =  6  case 
takes  approximately  10  minutes  of  1604  time. 

Accuracy:  Single  precision 

Caution  to  Users:  The  subroutine  as  written  requires 
that  the  data  series  should  contain  a  number  of  points 
exactly  a  power  of  two. 

Equipment  Configuration:  Standard  COOP 

References:  Writeup  of  subroutine  '!■  -■  G612COOL,  6/1/66 
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Writeup  of  program  TIES  Z24  SUBSET 


where 


T-l 

F^  (k)  =  x^t)  W  ~2 
t=0 

This  is  recognized  as  the  Fourier  transform  of  the  input 
data  computed  over  twice  its  length  with  zeros  filled  into 
the  second  half.  The  Cooley-Tukey  hyper-rapid  Fourier 
transform  routine  COOL  is  used  to  provide  the  high  speed 
necessary  here . 

Each  spectral  matrix  element  is  originally  T  +  1 
complex  points  long  between  DC  and  the  folding  frequency. 
It  is  then  smoothed  with  a  hanning  window  NS  times  to  its 
final  length  of  LF  points. 
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